{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mEnvironment variable CMDSTAN_HOME not found. Use set_cmdstan_home!.\u001b[39m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Turing]: AD chunk size is set as 40\n"
     ]
    }
   ],
   "source": [
    "using Turing, Distributions\n",
    "using PyPlot, PyCall\n",
    "using Mamba: describe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = [ 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 1.0, 1.0 ];\n",
    "N = length(y);  K = 3;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGgCAYAAABi2ofUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X100/XdP/7nJ2mb3qWBAm2SUqDclraATrQC6vCyUjkcZnf9psLXWXC4ne2UzYq39RLUeVN1c1M3LphOrV4OUS8Hbqggq7SOyc0Ae2lbbqVQKE3LbdKkbdomn98f5ZMSaaFpk7zzSZ6Pc3KOTT8Jr1RIn3m/X+/3W5JlWQYRERGRSmhEF0BERETkC4YXIiIiUhWGFyIiIlIVhhciIiJSFYYXIiIiUhWGFyIiIlIVhhciIiJSFYYXIiIiUhWGFyIiIlIVhhciIiJSFYYXIiIiUpUo0QX4g9vtxokTJ6DX6yFJkuhyiIiIqB9kWUZLSwvMZjM0mv6Pp4RFeDlx4gTS09NFl0FEREQDcOzYMYwcObLf14dFeNHr9QC6X3xSUpLgaoiIiKg/bDYb0tPTPb/H+ysswosyVZSUlMTwQkREpDK+tnywYZeIiIhUheGFiIiIVIXhhYiIiFSF4YWIiIhUheGFiIiIVIXhhYiIiFSF4YWIiIhUheGFiIiIVIXhhYiIiFSF4YWIiIhUxafwsmrVKkydOtWzDf+MGTPw6aefXvIxH3zwATIzMxEbG4spU6bgk08+8fq+LMtYsWIFTCYT4uLikJeXh4MHD/r+SoiIiCgi+BReRo4cieeeew67d+/Grl278B//8R+49dZbUVNT0+v1X375JRYuXIglS5bgq6++QkFBAQoKClBdXe255oUXXsArr7yC1atXY8eOHUhISEB+fj7a29sH98qIiIgoLEmyLMuDeYLk5GT85je/wZIlSy763h133AGHw4ENGzZ47rv22mtxxRVXYPXq1ZBlGWazGffffz8eeOABAIDVakVqairKysqwYMGCXv9Mp9MJp9Pp+Vo5ldJqtfJgRqIw02Rrx1tfHkFbp0t0KULdlJmK6yYMF10GkV/ZbDYYDAaff38P+FRpl8uFDz74AA6HAzNmzOj1mm3btmHZsmVe9+Xn52P9+vUAgLq6OlgsFuTl5Xm+bzAYkJubi23btvUZXkpLS/Hkk08OtHQiUpE/fn4I/7P9qOgyhPtw93FUrZgDjca303eJwpHP4eWbb77BjBkz0N7ejsTERKxbtw5ZWVm9XmuxWJCamup1X2pqKiwWi+f7yn19XdObkpISr1CkjLwQUfj5+vg5AMC8qSaMGRYvuBoxXvtnHWztXTh6phUZwxNEl0MknM/hZdKkSaiqqoLVasX//u//YtGiRaisrOwzwASCTqeDTqcL2p9HRGJ0udzYZ2kBADwwZ1LE/uLeevAU/u+4FTUnrBH7MyC6kM9LpWNiYjB+/HhcddVVKC0txbRp0/Dyyy/3eq3RaERTU5PXfU1NTTAajZ7vK/f1dQ0RRa5vTzrg7HIjUReF0cmROeoCAFlmAwCg5oRNcCVEoWHQ+7y43W6v5tkLzZgxA+Xl5V73bd682dMjk5GRAaPR6HWNzWbDjh07+uyjIaLIUXPCCgCYbNJHdK9Htrm7kZHhhaibT9NGJSUlmDt3LkaNGoWWlhasWbMGFRUV2LRpEwCgsLAQaWlpKC0tBQDce++9+P73v48XX3wR8+bNw9q1a7Fr1y68+uqrAABJklBcXIynn34aEyZMQEZGBpYvXw6z2YyCggI/v1QiUhvll3X2+ZGHSKWEl9oTVsiyDEmK3CBHBPgYXpqbm1FYWIjGxkYYDAZMnToVmzZtws033wwAqK+vh0bTM5gzc+ZMrFmzBo899hgeffRRTJgwAevXr0dOTo7nmoceeggOhwM/+9nPcO7cOVx33XXYuHEjYmNj/fQSiUitlJGXLHNkb4GQaUyCRgJO2TvQ3OJEahLfHymyDXqfl1Aw0HXiRBS6ZFnGtCc/g629C5/86vqIDzBzfl+JA012vLF4Ov4jM/XyDyBSgYH+/ubZRkQUko6fbYOtvQsxWg0mpCaKLkc4ZeqspoF9L0QML0QUkpQpo4nGRERr+VbFpl2iHnxHIKKQ5GnWNUV2s65CmTarabQKroRIPIYXIgpJ1Q3dv6Sz0yK710WhhLhjZ9pgbe0UXA2RWAwvRBSSepZJM7wAgCE+GiOHxgHg6AsRwwsRhZyTLU40tzghSd3LhKlbz34v7HuhyMbwQkQhR2nWzRiegASdz0ewha1sHhNABIDhhYhCEHfW7V3PiiNOG1FkY3ghopBTy36XXilh7tuTDrR3ugRXQyQOwwsRhRxlZIHhxVtqkg7DEmLgcsvYZ2kRXQ6RMAwvRBRSWto7ceR0KwBOG32XJEk9+71w6ogiGMMLEYWUvY3dIwomQyySE2IEVxN6ctLYtEvE8EJEIYVTRpfGYwKIGF6IKMRwpdGlKT+XfY02dLncgqshEoPhhYhCCnfWvbTRyfFI1EXB2eXG4VMO0eUQCcHwQkQhw9nlwsGm7p6X7DSOvPRGo5Ew2aQH0HP+E1GkYXghopBxwGJHl1vGkPhomA2xossJWdxplyIdwwsRhYwLm3UlSRJcTejicmmKdAwvRBQy2KzbPxce0CjLsuBqiIKP4YWIQgaXSffPhBQ9orUSbO1dOH62TXQ5REHH8EJEIcHllj0b1DG8XFpMlAYTU7ubdjl1RJGI4YWIQkLdKQfaOl2Ii9YiY3ii6HJCHjero0jG8EJEIUEZQcg06aHVsFn3crjiiCIZwwsRhYRabk7nk2yuOKIIxvBCRCGBK418M9mUBEkCmmxOnLI7RZdDFFQML0QknCzLXGnkowRdFDKGJwDg1BFFHoYXIhKu0dqOs62diNJInlU0dHk9fS+cOqLIwvBCRMIpIwfjUxIRG60VXI16cMURRSqGFyISTjlgkP0uvvGEFx7QSBGG4YWIhKvhSqMBUcLekdOtaGnvFFwNUfAwvBCRcLVs1h2Q5IQYmM6fvq3sTkwUCXwKL6Wlpbj66quh1+uRkpKCgoIC7N+//5KPmT17NiRJuug2b948zzWLFy++6Pu33HLLwF4REanKWUcHTljbAfSclkz9x/1eKBL5FF4qKytRVFSE7du3Y/Pmzejs7MScOXPgcDj6fMxf//pXNDY2em7V1dXQarW47bbbvK675ZZbvK579913B/aKiEhVlCmj0cPioY+NFlyN+mRxp12KQFG+XLxx40avr8vKypCSkoLdu3fjhhtu6PUxycnJXl+vXbsW8fHxF4UXnU4Ho9HYrzqcTieczp5NmWw2/qMlUivu7zI4XHFEkWhQPS9Wa/ebzncDyqW8/vrrWLBgARISErzur6ioQEpKCiZNmoRf/OIXOH36dJ/PUVpaCoPB4Lmlp6cP7AUQkXDcWXdwlPBysKkFzi6X4GqIgmPA4cXtdqO4uBizZs1CTk5Ovx6zc+dOVFdX45577vG6/5ZbbsHbb7+N8vJyPP/886isrMTcuXPhcvX+D7GkpARWq9VzO3bs2EBfBhEJpoy8sN9lYNKGxMEQF40ut4yDTXbR5RAFhU/TRhcqKipCdXU1tm7d2u/HvP7665gyZQquueYar/sXLFjg+e8pU6Zg6tSpGDduHCoqKnDTTTdd9Dw6nQ46nW6gpRNRiGjt6MLhU909c5w2GhhJkpCTloR/HTqNmhNW5KRxBIvC34BGXpYuXYoNGzZgy5YtGDlyZL8e43A4sHbtWixZsuSy144dOxbDhw/HoUOHBlIeEanE3sYWyDIwQq9Dij5WdDmqlc2mXYowPo28yLKMX/7yl1i3bh0qKiqQkZHR78d+8MEHcDqd+PGPf3zZa48fP47Tp0/DZDL5Uh4RqQz3d/EPNu1SpPFp5KWoqAjvvPMO1qxZA71eD4vFAovFgra2Ns81hYWFKCkpueixr7/+OgoKCjBs2DCv++12Ox588EFs374dR44cQXl5OW699VaMHz8e+fn5A3xZRKQG3FnXP5Sf395GG1xuWXA1RIHn08jLqlWrAHRvPHehN998E4sXLwYA1NfXQ6PxzkT79+/H1q1b8dlnn130nFqtFl9//TXeeustnDt3DmazGXPmzMFTTz3FvhaiMFd9fuQlhyuNBiVjeCLiorVo7XCh7pQD41MSRZdEFFA+TxtdTkVFxUX3TZo0qc/HxsXFYdOmTb6UQURhoNPlxgFL9+oYLpMeHK1GQqZJj6/qz6HmhJXhhcIezzYiIiEONtnR4XJDHxuF9OQ40eWonjJ1VMu+F4oADC9EJIRnfxdTEiRJElyN+nHFEUUShhciEoI76/rXhQc09meKn0jNGF6ISIharjTyq4mpemg1Es62dqLx/CndROGK4YWIgs7tllHbeD68pDG8+ENstBYTzjfqcuqIwh3DCxEFXf2ZVtidXYiJ0mDcCK6M8ZesC6aOiMIZwwsRBZ0yMpBp1CNay7chf2HTLkUKvmsQUdDV8FiAgMjhcmmKEAwvRBR0yshAFlca+ZUybdRwrg1nHR2CqyEKHIYXIgo6nmkUGPrYaIweFg8AnoZoonDE8EJEQdVsa8cpuxMaCZhsZHjxt2w27VIEYHghoqBSDmMcNyIRcTFawdWEH6Vpt7qBIy8UvhheiCioaho4ZRRIXC5NkYDhhYiCiscCBJYSCg+fcqC1o0twNUSBwfBCREFV08hl0oGUoo/FCL0OsgzsbWwRXQ5RQDC8EFHQWNs6cexMG4Ce6Q3yv2zPfi+cOqLwxPBCREGjbJ6WNiQOQ+JjBFcTvnpWHLFpl8ITwwsRBQ131g0OHhNA4Y7hhYiCppbNukGhhMP9lhZ0utyCqyHyP4YXIgoa7qwbHKOS46GPjUKHy41DzXbR5RD5HcMLEQVFe6cLh052/yLNTmN4CSRJkpBlYt8LhS+GFyIKiv2WFrjcMpITYmBMihVdTtjr6XvhiiMKPwwvRBQUF04ZSZIkuJrwxxVHFM4YXogoKJQzjbi/S3AoU3O1J2xwu2XB1RD5F8MLEQUFjwUIrnEjEhETpYHd2YX6M62iyyHyK4YXIgq4Lpcb+xq7w0sOR16CIlqrQaZRD4BTRxR+GF6IKOAOn3LA2eVGQowWY4YliC4nYmTzhGkKUwwvRBRwyi/PyaYkaDRs1g2WLO60S2GK4YWIAq6mgZvTicAVRxSuGF6IKODYrCvGZGMSNBJwyu5Es61ddDlEfsPwQkQBJcuyZ9qIy6SDKy5Gi7EjEgFw9IXCi0/hpbS0FFdffTX0ej1SUlJQUFCA/fv3X/IxZWVlkCTJ6xYb6727pizLWLFiBUwmE+Li4pCXl4eDBw/6/mqIKOQcP9sGW3sXorUSJqbqRZcTcdi0S+HIp/BSWVmJoqIibN++HZs3b0ZnZyfmzJkDh8NxycclJSWhsbHRczt69KjX91944QW88sorWL16NXbs2IGEhATk5+ejvZ3DnERqp3zin5CiR0wUB3uDLYdNuxSGony5eOPGjV5fl5WVISUlBbt378YNN9zQ5+MkSYLRaOz1e7Is46WXXsJjjz2GW2+9FQDw9ttvIzU1FevXr8eCBQsueozT6YTT6fR8bbPxHyVRqKo9/4mfzbpisGmXwtGgPgZZrd1vSsnJyZe8zm63Y/To0UhPT8ett96Kmpoaz/fq6upgsViQl5fnuc9gMCA3Nxfbtm3r9flKS0thMBg8t/T09MG8DCIKoAvPNKLgU/qM6s+0wtbeKbgaIv8YcHhxu90oLi7GrFmzkJOT0+d1kyZNwhtvvIGPPvoI77zzDtxuN2bOnInjx48DACwWCwAgNTXV63Gpqame731XSUkJrFar53bs2LGBvgwiCjBPeEnjSiMRhsTHIG1IHIDuc46IwoFP00YXKioqQnV1NbZu3XrJ62bMmIEZM2Z4vp45cyYmT56MP/3pT3jqqacG9GfrdDrodLoBPZaIgueU3QmLrR2S1L1BHYmRbU5Cw7k2VDdYce3YYaLLIRq0AY28LF26FBs2bMCWLVswcuRInx4bHR2NK6+8EocOHQIATy9MU1OT13VNTU199skQkToooy5jhiUgUTfgz0o0SMr+Ohx5oXDhU3iRZRlLly7FunXr8PnnnyMjI8PnP9DlcuGbb76ByWQCAGRkZMBoNKK8vNxzjc1mw44dO7xGbIhIfWrYrBsS2LRL4canj0JFRUVYs2YNPvroI+j1ek9PisFgQFxc95xqYWEh0tLSUFpaCgD49a9/jWuvvRbjx4/HuXPn8Jvf/AZHjx7FPffcA6B7JVJxcTGefvppTJgwARkZGVi+fDnMZjMKCgr8+VqJKMi4s25oyE7rDi+HTtrR3ulCbLRWcEVEg+NTeFm1ahUAYPbs2V73v/nmm1i8eDEAoL6+HhpNz4DO2bNn8dOf/hQWiwVDhw7FVVddhS+//BJZWVmeax566CE4HA787Gc/w7lz53Dddddh48aNF21mR0TqUsuVRiHBmBSL5IQYnHF0YL+lBdPSh4guiWhQJFmWZdFFDJbNZoPBYIDVakVSEt8kiUKB3dmFnMc3AQB2P5aHYYlsshfprtd34J8HT+HZH07B/8sdJbocIgAD//3N7S6JKCD2NnaPuhiTYhlcQkAWjwmgMMLwQkQBUdPAZt1Qks1jAiiMMLwQUUBwZ93QknP+/8M+iw0ut+q7BSjCMbwQUUAo4SWLK41CwphhCUiI0aK9043DJ+2iyyEaFIYXIvK7ji43Dja3AODIS6jQaCTPLsecOiK1Y3ghIr870NSCTpcMQ1w0Rg6NE10OnZfNpl0KEwwvROR3yi/HLFMSJEkSXA0plKbd6gaOvJC6MbwQkd+xWTc0XbhcOgy2+KIIxvBCRH7nCS9pDC+hZGKqHtFaCbb2Lhw/2ya6HKIBY3ghIr9yuWXPBnU80yi0xERpMCFFD4BNu6RuDC9E5FdHTjvQ2uFCbLQGY4cniC6HvkOZyqtl0y6pGMMLEfmV8ok+05iEKC3fYkJNz4ojjryQevGdhYj8SllpxGbd0JSdxmMCSP0YXojIr2pPsN8llE02JUGSAIutHaftTtHlEA0IwwsR+Y0sy1wmHeISdVEYM6y7F4mjL6RWDC9E5DcWWzvOODqg1UiYZNSLLof6wL4XUjuGFyLym5rzO7eOH5GI2Git4GqoL8qUHo8JILVieCEiv+GUkTr0LJfmyAupE8MLEfmN50wjhpeQpoSXutMO2J1dgqsh8h3DCxH5TQ1XGqnCsEQdjEmxkGV4dkMmUhOGFyLyi7OODjSc6z4vhyMvoc/TtNvAvhdSH4YXIvKL2vOf4NOT42CIixZcDV0OVxyRmjG8EJFfeHbWNXHKSA2yzNxpl9SL4YWI/IIrjdRF+f90sLkFHV1uwdUQ+YbhhYj8QgkvOWkceVGDkUO7p/c6XTIONLWILofIJwwvRDRobR0uHD5pB8CRF7WQJAlZJu73QurE8EJEg7bXYoNbBoYn6pCSFCu6HOqnnqZdrjgidWF4IaJBY7+LOilTfGzaJbVheCGiQatVVhoxvKiK8v9rb6MNbrcsuBqi/mN4IaJB48666jR2RCJiozVwdLhw5LRDdDlE/cbwQkSD0ulyY5+le7UKR17URauRkGnkZnWkPj6Fl9LSUlx99dXQ6/VISUlBQUEB9u/ff8nHvPbaa7j++usxdOhQDB06FHl5edi5c6fXNYsXL4YkSV63W265xfdXQ0RBd6jZjo4uNxJ1URiVHC+6HPKREjir2bRLKuJTeKmsrERRURG2b9+OzZs3o7OzE3PmzIHD0fdwY0VFBRYuXIgtW7Zg27ZtSE9Px5w5c9DQ0OB13S233ILGxkbP7d133x3YKyKioFI+sWeZkqDRSIKrIV8pU31cLk1qEuXLxRs3bvT6uqysDCkpKdi9ezduuOGGXh/zl7/8xevrP//5z/jwww9RXl6OwsJCz/06nQ5Go7FfdTidTjidTs/XNhv/0RGJoiyz5WGM6nThGUeyLEOSGEAp9A2q58Vq7X7TSk5O7vdjWltb0dnZedFjKioqkJKSgkmTJuEXv/gFTp8+3edzlJaWwmAweG7p6ekDewFENGhcJq1uk4x6aDUSzjg6YLG1iy6HqF8GHF7cbjeKi4sxa9Ys5OTk9PtxDz/8MMxmM/Ly8jz33XLLLXj77bdRXl6O559/HpWVlZg7dy5cLlevz1FSUgKr1eq5HTt2bKAvg4gGwe2WsZcrjVQtNlqL8SMSAQA1DRzFJnXwadroQkVFRaiursbWrVv7/ZjnnnsOa9euRUVFBWJje3bhXLBggee/p0yZgqlTp2LcuHGoqKjATTfddNHz6HQ66HS6gZZORH5y7GwrWpxdiNFqMCE1UXQ5NEDZ5iTsb2pBzQkb8rJSRZdDdFkDGnlZunQpNmzYgC1btmDkyJH9esxvf/tbPPfcc/jss88wderUS147duxYDB8+HIcOHRpIeUQUJMqU0SSjHtFa7rygVlk8JoBUxqeRF1mW8ctf/hLr1q1DRUUFMjIy+vW4F154Ac888ww2bdqE6dOnX/b648eP4/Tp0zCZTL6UR0RBVsOddcOCMuXHvV5ILXz6qFRUVIR33nkHa9asgV6vh8VigcViQVtbm+eawsJClJSUeL5+/vnnsXz5crzxxhsYM2aM5zF2e/cJtHa7HQ8++CC2b9+OI0eOoLy8HLfeeivGjx+P/Px8P71MIgoENuuGB2XkpeFcG861dgiuhujyfAovq1atgtVqxezZs2EymTy39957z3NNfX09GhsbvR7T0dGBH/3oR16P+e1vfwsA0Gq1+Prrr/GDH/wAEydOxJIlS3DVVVfhn//8J/taiEKcZ48XNuuqmiEu2rPBIPd7ITXwedrocioqKry+PnLkyCWvj4uLw6ZNm3wpg4hCQHNLO062OCFJwGSTXnQ5NEjZ5iTUn2lFzQkbZo4fLrocoktihx0RDYgy6jJ2eALiYwa8cJFCRDabdklFGF6IaEBqub9LWGHTLqkJwwsRDUh1A1cahRPl/+O3J+1o6+h9g1CiUMHwQkQDUsORl7CSkhSL4Yk6uGVgr4WjLxTaGF6IyGe29k7Un2kFwJGXcHLhIY1EoYzhhYh8pvS7mA2xGJoQI7ga8hclvNSyaZdCHMMLEfmM+7uEJzbtklowvBCRz3gsQHhS/n/us7Sg0+UWXA1R3xheiMhntTwWICyNSo5Hoi4KHV1ufHvSLrocoj4xvBCRT9o7XTjY3P2LLSeN00bhRKORkGU637TbwKkjCl0ML0TkkwNNLXC5ZQyNj4bJECu6HPKzLK44IhVgeCEin1y4v4skSYKrIX9TRtN4TACFMoYXIvIJm3XDm2e5dKOtX4fxEonA8EJEPulZJs3wEo7GpyQiJkqDlvYuHDvTJrocol4xvBBRv7ncMvY28liAcBat1WBSqh4AUM2pIwpRDC9E1G+HT9rR3ulGXLQWGcMTRJdDAdJzTADDC4Umhhci6jdlymiySQ+ths264YpnHFGoY3ghon7radbllFE4y+IxARTiGF6IqN9quLNuRJhs0kOSgJMtTjS3tIsuh+giDC9E1C+yLHvt8ULhKz4mCmPP9zRx9IVCEcMLEfVLw7k2WNs6EaWRMNGYKLocCjAloNYyvFAIYnghon5RPoGPT0mELkoruBoKNK44olDG8EJE/aKEFx7GGBmy2bRLIYzhhYj6pZbHAkQU5f/z0dOtsLV3Cq6GyBvDCxH1C5t1I8vQhBikDYkDAOzl6AuFGIYXIrqsM44ONFq7l8xONukFV0PBksXN6ihEMbwQ0WUpTZtjhsVDHxstuBoKFu60S6GK4YWILqu6gVNGkainaZcrjii0MLwQ0WUpv7yy2KwbUZSRl4PNdrR3ugRXQ9SD4YWILquWxwJEJJMhFkPjo+FyyzjQ1CK6HCIPhhciuiSHswt1px0AOG0UaSRJ4n4vFJJ8Ci+lpaW4+uqrodfrkZKSgoKCAuzfv/+yj/vggw+QmZmJ2NhYTJkyBZ988onX92VZxooVK2AymRAXF4e8vDwcPHjQt1dCRAGxt9EGWQZS9DqM0OtEl0NBxp12KRT5FF4qKytRVFSE7du3Y/Pmzejs7MScOXPgcDj6fMyXX36JhQsXYsmSJfjqq69QUFCAgoICVFdXe6554YUX8Morr2D16tXYsWMHEhISkJ+fj/Z2nmZKJBpPko5sXC5NoUiSZVke6INPnjyJlJQUVFZW4oYbbuj1mjvuuAMOhwMbNmzw3HfttdfiiiuuwOrVqyHLMsxmM+6//3488MADAACr1YrU1FSUlZVhwYIFFz2n0+mE0+n0fG2z2ZCeng6r1YqkJL7Bkv/Isox3th/F4VN9B/Rwt+PwGdQ22rD0xvF4IH+S6HIoyA4125H3u0roojT4f7mjRJcjTFy0FotnjkFKUqzoUsKKzWaDwWDw+fd31GD+UKu1exgxOTm5z2u2bduGZcuWed2Xn5+P9evXAwDq6upgsViQl5fn+b7BYEBubi62bdvWa3gpLS3Fk08+OZjSifrl6+NWLP+oRnQZIWFa+hDRJZAAGcMTkBQbBVt7F9781xHR5QjV1unC4/OzRZdBGER4cbvdKC4uxqxZs5CTk9PndRaLBampqV73paamwmKxeL6v3NfXNd9VUlLiFYiUkRcif/v6+DkAwLgRCbglxyi4GnFS9LG4KTNFdBkkgFYj4bXC6fji4EnRpQhTd8qBT76x4Ovj7PsJFQMOL0VFRaiursbWrVv9WU+/6HQ66HRsHKTAU+b587ONeDA/U3A1RGLkjh2G3LHDRJchzIGmFnzyjQV7G21wuWVoNZLokiLegJZKL126FBs2bMCWLVswcuTIS15rNBrR1NTkdV9TUxOMRqPn+8p9fV1DJAoPIySiscMToIvSoLXDhSOnI7f/LZT4FF5kWcbSpUuxbt06fP7558jIyLjsY2bMmIHy8nKv+zZv3owZM2YAADIyMmA0Gr2usdls2LFjh+caIhE6XW7st3RvzMWVNkSRK0qrQaaJq65CiU/hpaioCO+88w7WrFkDvV4Pi8UCi8WCtrY2zzWFhYUoKSnxfH3vvfdi48aNePHFF7Fv3z488cQT2LVrF5YuXQqgexOk4uJiPP0xyOZ1AAAgAElEQVT00/jb3/6Gb775BoWFhTCbzSgoKPDTyyTy3aFmOzpcbiTqojAqOV50OUQkEPe7CS0+9bysWrUKADB79myv+998800sXrwYAFBfXw+NpicTzZw5E2vWrMFjjz2GRx99FBMmTMD69eu9mnwfeughOBwO/OxnP8O5c+dw3XXXYePGjYiN5ZI0Ekf5hJVlSoKGc9xEEU0JL7UceQkJPoWX/mwJU1FRcdF9t912G2677bY+HyNJEn7961/j17/+tS/lEAWU8gkrO41TRkSRLueCYxJkWYYk8QONSDzbiKgPbNYlIsUkox5ajYQzjg5YbNz9XTSGF6JeuN0y9nJbfCI6LzZai/EjEgEANQ2cOhKN4YWoF8fOtqLF2YWYKA3GpySKLoeIQkA2z3kKGQwvRL2oPv/JalKqHtFa/jMhop5DKqu54kg4visT9cLTrMspIyI6T+l/44oj8RheiHpRw34XIvoOZeSl4Vwbzjo6BFcT2RheiHrh2eOFK42I6DxDXDTSk+MAALWNHH0RieGF6Duabe04ZXdCkoDJJr3ocogohGSblP1e2PciEsML0Xcooy5jhycgPmbAB68TURjiiqPQwPBC9B09zbqcMiIib8qO2wwvYjG8EH0Hm3WJqC/Kh5rDJ+1o63AJriZyMbwQfYcSXnLSOPJCRN5S9DoMT9TBLQN7LRx9EYXhhegCtvZO1J9pBcCRFyK6mCRJ7HsJAQwvRBdQNp9KGxKHIfExgqsholCkhJdarjgShuGF6AI9+7tw1IWIeqf0vXDkRRyGF6IL8FgAIroc5f1hn6UFnS634GoiE8ML0QVqPSuN2KxLRL0blRyPRF0UOrrc+PakXXQ5EYnhhei89k4XDjZ3vxFx5IWI+qLRSMgynT9huoFTRyIwvBCdt9/SApdbxtD4aJgMsaLLIaIQluVZccSmXREYXojOq7lgykiSJMHVEFEo43JpsRheiM5jsy4R9ZfSF7f3hA1utyy4msjD8EJ0HpdJE1F/TUhNRIxWgxZnF46dbRVdTsRheCEC4HLL2GfhSiMi6p9orQYTjYkAOHUkAsMLEboPWWvvdCMuWouM4QmiyyEiFcg2KZvVsWk32BheiOA9ZaTVsFmXiC4vJ41Nu6IwvBCBzbpE5LssHhMgDMMLES5cJs3wQkT9M9mkhyQBJ1ucaG5pF11ORGF4oYgny7LXHi9ERP0RHxOFsed75Dj6ElwMLxTxGs61wdrWiSiNhAmpiaLLISIVUT7w1DK8BBXDC0U85RPThFQ9dFFawdUQkZpk85gAIRheKOLVNLBZl4gGRhl54QGNweVzePniiy8wf/58mM1mSJKE9evXX/L6xYsXQ5Kki27Z2dmea5544omLvp+Zmen7qyEaADbrEtFAKe8b9WdaYWvvFFxN5PA5vDgcDkybNg0rV67s1/Uvv/wyGhsbPbdjx44hOTkZt912m9d12dnZXtdt3brV19KIBoTNukQ0UEMTYmA+fwo9+16CJ8rXB8ydOxdz587t9/UGgwEGQ88vhfXr1+Ps2bO4++67vQuJioLRaOzXczqdTjidTs/XNhv/wtDAnLY7YbF1L3GcbNILroaI1CjLbMAJaztqTthw7dhhosuJCEHveXn99deRl5eH0aNHe91/8OBBmM1mjB07FnfeeSfq6+v7fI7S0lJPKDIYDEhPTw902RSmlFGXMcPioY+NFlwNEakRm3aDL6jh5cSJE/j0009xzz33eN2fm5uLsrIybNy4EatWrUJdXR2uv/56tLS09Po8JSUlsFqtntuxY8eCUT6FIU4ZEdFgKeGF00bB4/O00WC89dZbGDJkCAoKCrzuv3AaaurUqcjNzcXo0aPx/vvvY8mSJRc9j06ng06nC3i9FP48xwKksVmXiAYmJ637w8/BZjvaO12IjeaWC4EWtJEXWZbxxhtv4K677kJMTMwlrx0yZAgmTpyIQ4cOBak6ilS1HHkhokEyGWIxND4aLreMA029zxiQfwUtvFRWVuLQoUO9jqR8l91ux7fffguTyRSEyihSOZxdqDvtAMBl0kQ0cJIkeT4A8ZiA4PA5vNjtdlRVVaGqqgoAUFdXh6qqKk+DbUlJCQoLCy963Ouvv47c3Fzk5ORc9L0HHngAlZWVOHLkCL788kv88Ic/hFarxcKFC30tj6jf9jbaIMtAapIOwxM5DUlEA8em3eDyuedl165duPHGGz1fL1u2DACwaNEilJWVobGx8aKVQlarFR9++CFefvnlXp/z+PHjWLhwIU6fPo0RI0bguuuuw/bt2zFixAhfyyPqNzbrEpG/ZHnCC0degsHn8DJ79mzIstzn98vKyi66z2AwoLW1tc/HrF271tcyiAbN06zLKSMiGiTlQ9C+xha43DK0GklwReGNZxtRxOKxAETkLxnDExAXrUVbpwt1p+yiywl7DC8UkTq63J5VAZw2IqLB0mokzy7dnDoKPIYXikgHmlrQ6ZKRFBuFkUPjRJdDRGGg54RpNu0GGsMLRSRlf5cscxIkiXPTRDR42WzaDRqGF4pIPc26nDIiIv+4cK+XSy1socFjeKGIxGZdIvK3icZERGkkWNs60XCuTXQ5YY3hhSKO2y1jbyP3eCEi/9JFaTE+JREAp44CjeGFIs6R0w44OlzQRWkwbkSC6HKIKIwohzQyvAQWwwtFHOVNJdOUhCgt/wkQkf8oU9G1PCYgoPjOTRGH/S5EFCg8oDE4GF4o4vBYACIKFGWjukZrO844OgRXE74YXiiiyLLs2eOFzbpE5G/62GiMGRYPgCdMBxLDC0WUJpsTpx0d0GokZBr1osshojDEqaPAY3ihiKJ8Eho3IgGx0VrB1RBROMriTrsBx/BCEaWGU0ZEFGA9xwRw2ihQGF4oorBZl4gCTflwVHfKAYezS3A14YnhhSJKdUPPgYxERIEwQq9Dil4HWYZnN2/yL4YXihjnWjs8541kmzhtRESBwxOmA4vhhSKGskR65NA4GOKjBVdDROGsZ8UR+14CgeGFIobyCSiHzbpEFGA5aRx5CSSGF4oYbNYlomBRRl4ONLWgo8stuJrww/BCEcOzTDqN4YWIAmvk0DgkxUah0yXjYHOL6HLCDsMLRYS2Dhe+PWkHwD1eiCjwJEniZnUBxPBCEWGfxQa3DAxPjEGKXie6HCKKAMoHpVqGF79jeKGIoHzyyTIbIEmS4GqIKBJwp93AYXihiNBzLAD7XYgoOC4ceXG7ZcHVhBeGF4oItVxpRERBNm5EAnRRGjg6XDh6plV0OWGF4YXCXpfLjX2W7m5/NusSUbBEaTXINOoBcOrI3xheKOx9e9IBZ5cbiboojE6OF10OEUWQLM9Ou2za9SeGFwp71Q3dn3gmm/TQaNisS0TBo0xVK+9D5B8MLxT2epp1OWVERMGlhJfaEzbIMpt2/cXn8PLFF19g/vz5MJvNkCQJ69evv+T1FRUVkCTpopvFYvG6buXKlRgzZgxiY2ORm5uLnTt3+loaUa+UueYsNusSUZBlGpOgkYDTjg402ZyiywkbPocXh8OBadOmYeXKlT49bv/+/WhsbPTcUlJSPN977733sGzZMjz++OPYs2cPpk2bhvz8fDQ3N/taHpEXWZZR28gDGYlIjLgYLcanJAJg064/Rfn6gLlz52Lu3Lk+/0EpKSkYMmRIr9/73e9+h5/+9Ke4++67AQCrV6/Gxx9/jDfeeAOPPPLIRdc7nU44nT0J1mZjIxT17tiZNrS0dyFGq8GE1ETR5RBRBMo2G3CgyY6aEzbcNDlVdDlhIWg9L1dccQVMJhNuvvlm/Otf//Lc39HRgd27dyMvL6+nKI0GeXl52LZtW6/PVVpaCoPB4Lmlp6cHvH5SJ+WTzkRjIqK1bPEiouDjTrv+F/B3c5PJhNWrV+PDDz/Ehx9+iPT0dMyePRt79uwBAJw6dQoulwupqd5pNDU19aK+GEVJSQmsVqvnduzYsUC/DFIpT7OuiVNGRCQGD2j0P5+njXw1adIkTJo0yfP1zJkz8e233+L3v/89/ud//mdAz6nT6aDT8XA9ujzlk052Gpt1iUgM5cPT8bNtsLZ2whAfLbgi9RMyjn7NNdfg0KFDAIDhw4dDq9WiqanJ65qmpiYYjUYR5VEY4ZlGRCSaIT4aI4fGAQBqGjl15A9CwktVVRVMJhMAICYmBldddRXKy8s933e73SgvL8eMGTNElEdh4mSLE80tTkhS93JFIiJRLtzvhQbP52kju93uGTUBgLq6OlRVVSE5ORmjRo1CSUkJGhoa8PbbbwMAXnrpJWRkZCA7Oxvt7e3485//jM8//xyfffaZ5zmWLVuGRYsWYfr06bjmmmvw0ksvweFweFYfEQ2EMmWUMTwBCbqAz5ASEfUp22zAppom9r34ic/v6Lt27cKNN97o+XrZsmUAgEWLFqGsrAyNjY2or6/3fL+jowP3338/GhoaEB8fj6lTp+If//iH13PccccdOHnyJFasWAGLxYIrrrgCGzduvKiJl8gX3FmXiEIFVxz5lySHwX7FNpsNBoMBVqsVSUmcHqBuRX/Zg4+/acQjczPx8++PE10OEUUwi7Ud15aWQ6uRUPNkPmKjtaJLCgkD/f3NjS8obHlWGrFZl4gES03SYVhCDFxuGfssLaLLUT2GFwpLtvZOHDndCoDTRkQkniRJnv1eeML04DG8UFjae77fxWyIRXJCjOBqiIiAnLTuD1Js2h08hhcKS8qbQxZHXYgoRPQsl+bIy2AxvFBY4uZ0RBRqlCnsfZYWdLncgqtRN4YXCkts1iWiUDM6OR6Juig4u9z49qRDdDmqxvBCYcfZ5cKhZjsAIDuN00ZEFBo0GgmTTXoA3O9lsBheKOwcsNjR5ZYxJD4aZkOs6HKIiDyUqSM27Q4OwwuFnQunjCRJElwNEVGPLO606xcMLxR2eCwAEYWqCw9oDIMN7oVheKGww2ZdIgpVE1L0iNZKsLV34fjZNtHlqBbDC4UVl1vG3sburbcZXogo1MREaTAxlU27g8XwQmGl7pQDbZ0uxEVrkTE8UXQ5REQX6Tlhmk27A8XwQmFF+SSTadJDq2GzLhGFHq44GjyGFwor3FmXiEJdNg9oHDSGFworyshLDlcaEVGImmxKgiQBzS1OnGxxii5HlRheKGzIssxl0kQU8hJ0UcgYngCATbsDxfBCYeOEtR3nWjsRpZEw0chmXSIKXex7GRyGFwobNefnj8enJEIXpRVcDRFR3y7crI58x/BCYYNTRkSkFtk8JmBQGF4obHClERGphfIh68jpVrS0dwquRn0YXihs1PJYACJSieSEGJjOn3qv7ApO/cfwQmHhrKMDJ6ztAHpObSUiCmWcOho4hhcKC8qU0ehh8dDHRguuhojo8rK44mjAGF4oLPAkaSJSG55xNHAMLxQWuNKIiNRGCS8Hm1rg7HIJrkZdGF4oLCgjL+x3ISK1SBsSB0NcNLrcMg422UWXoyoML6R6rR1dOHzKAYBnGhGRekiShJw0Nu0OBMMLqd7eRhtkGUjR6zBCrxNdDhFRvylT3dUN7HvxBcMLqR43pyMiteJy6YFheCHVq2lgsy4RqZMSXvY2tsDllgVXox4+h5cvvvgC8+fPh9lshiRJWL9+/SWv/+tf/4qbb74ZI0aMQFJSEmbMmIFNmzZ5XfPEE09AkiSvW2Zmpq+lUYSqaeQyaSJSp4zhiYiL1qKt04W68717dHk+hxeHw4Fp06Zh5cqV/br+iy++wM0334xPPvkEu3fvxo033oj58+fjq6++8rouOzsbjY2NntvWrVt9LY0iUKfLjQOW7i59jrwQkdpoNRIyTXoAnDryRZSvD5g7dy7mzp3b7+tfeuklr6+fffZZfPTRR/j73/+OK6+8sqeQqCgYjcZ+PafT6YTT6fR8bbOx0SlSHWyyo8Plhj42CunJcaLLISLyWbY5CV/Vn0PtCRtuvSJNdDmqEPSeF7fbjZaWFiQnJ3vdf/DgQZjNZowdOxZ33nkn6uvr+3yO0tJSGAwGzy09PT3QZVOI8uzvYkqCJEmCqyEi8l02jwnwWdDDy29/+1vY7Xbcfvvtnvtyc3NRVlaGjRs3YtWqVairq8P111+PlpbeT9osKSmB1Wr13I4dOxas8inEcGddIlK7C1ccyTKbdvvD52mjwVizZg2efPJJfPTRR0hJSfHcf+E01NSpU5Gbm4vRo0fj/fffx5IlSy56Hp1OB52O+3kQUMtl0kSkchNT9dBqJJxt7USjtR3mIZwCv5ygjbysXbsW99xzD95//33k5eVd8tohQ4Zg4sSJOHToUJCqIzVyu2XUNp4PL2kML0SkTrHRWkxISQTAqaP+Ckp4effdd3H33Xfj3Xffxbx58y57vd1ux7fffguTyRSE6kit6s+0wu7sQkyUBuNGJIouh4howLK4WZ1PfA4vdrsdVVVVqKqqAgDU1dWhqqrK02BbUlKCwsJCz/Vr1qxBYWEhXnzxReTm5sJiscBiscBq7fkf9MADD6CyshJHjhzBl19+iR/+8IfQarVYuHDhYF8fhTHlE0qmUY9oLfdbJCL1YtOub3x+x9+1axeuvPJKzzLnZcuW4corr8SKFSsAAI2NjV4rhV599VV0dXWhqKgIJpPJc7v33ns91xw/fhwLFy7EpEmTcPvtt2PYsGHYvn07RowYMdjXR2FM+YTCZl0iUruc8yMvtQwv/eJzw+7s2bMv2Q1dVlbm9XVFRcVln3Pt2rW+lkGEajbrElGYUKaNGs614ayjA0MTYgRXFNo41k6qJMsyak/wWAAiCg/62GiMHhYPgFNH/cHwQqrU3OLEKXsHNBKQaWR4ISL14wnT/cfwQqqk/OMeNyIRcTFawdUQEQ0em3b7j+GFVKmmgf0uRBReuFy6/xheSJV4LAARhRvlw9jhUw60dnQJria0MbyQKtU0slmXiMJLij4WI/Q6yDKwt7H3s/2oG8MLqY61rRPHzrQB6BlmJSIKB9me/V44dXQpDC+kOsomTmlD4jAknnshEFH46FlxxKbdS2F4IdWp4f4uRBSmuOKofxheSHVq2axLRGFK+VC239KCTpdbcDWhi+GFVKeGxwIQUZhKHxoPvS4KHS43DjXbRZcTshheSFXaO104dLL7H3ROGkdeiCi8aDTSBfu9cOqoLwwvpCr7LS1wuWUMS4hBapJOdDlERH7X0/fCFUd9YXghVak+/485y5wESZIEV0NE5H+eFUcNHHnpC8MLqQp31iWicJeddn6vl0Yb3G5ZcDWhieGFVIXNukQU7saNSERMlAZ2Zxfqz7SKLickMbyQanS53NjXyPBCROEtWqtBplEPgE27fWF4IdU4fMoBZ5cbCTFajBmWILocIqKAyeYJ05fE8EKqofwjnmxKgkbDZl0iCl9Z3Gn3khheSDWUzntOGRFRuOMZR5fG8EKqwZVGRBQpJhuToJGAU3Ynmm3tossJOQwvpAqyLHumjbI48kJEYS4uRouxIxIBcPSlNwwvpArHz7bB1t6FaK2Eial60eUQEQUcm3b7xvBCqqB88piYqkdMFP/aElH4y2HTbp/4W4BUofb8Jw826xJRpGDTbt8YXkgV2KxLRJFG6e+rP9MKa1un4GpCC8MLqUI1R16IKMIMiY9B2pA4AEAtR1+8MLxQyDtld6LJ5oQkdW9QR0QUKdi02zuGFwp5ypRRxrAEJOiiBFdDRBQ8ylQ5R168MbxQyOP+LkQUqdi02zuGFwp5bNYlokiVndYdXg6dtKO90yW4mtDhc3j54osvMH/+fJjNZkiShPXr11/2MRUVFfje974HnU6H8ePHo6ys7KJrVq5ciTFjxiA2Nha5ubnYuXOnr6VRmKo9wTONiCgyGZNikZwQA5dbxn5Li+hyQobP4cXhcGDatGlYuXJlv66vq6vDvHnzcOONN6KqqgrFxcW45557sGnTJs817733HpYtW4bHH38ce/bswbRp05Cfn4/m5mZfy6MwY3d2oe6UAwDDCxFFHkmSOHXUC0mWZXnAD5YkrFu3DgUFBX1e8/DDD+Pjjz9GdXW1574FCxbg3Llz2LhxIwAgNzcXV199Nf74xz8CANxuN9LT0/HLX/4SjzzyyEXP6XQ64XQ6PV/bbDakp6fDarUiKcl/v+C6XG4888levz0f+e6sowPrq07AmBSL7Y/eJLocIqKgK/10L/5UeRhT0gyYPmaosDqiNBL+a16WX5/TZrPBYDD4/Ps74Es3tm3bhry8PK/78vPzUVxcDADo6OjA7t27UVJS4vm+RqNBXl4etm3b1utzlpaW4sknnwxc0ee5ZeDNfx0J+J9Dlzctnf0uRBSZrkwfAgD4psGKbxrELZmOidL4PbwMVMDDi8ViQWpqqtd9qampsNlsaGtrw9mzZ+FyuXq9Zt++fb0+Z0lJCZYtW+b5Whl58TeNBBTdOM7vz0u+idZq8P99b6ToMoiIhLg5y4gn5mfhpN15+YsDSKsJnTU+qtw0Q6fTQafTBfzPidJq8GB+ZsD/HCIior5oNRIWz8oQXUZICXh4MRqNaGpq8rqvqakJSUlJiIuLg1arhVar7fUao9EY6PKIiIhIZQI+BjRjxgyUl5d73bd582bMmDEDABATE4OrrrrK6xq3243y8nLPNUREREQKn8OL3W5HVVUVqqqqAHQvha6qqkJ9fT2A7n6UwsJCz/U///nPcfjwYTz00EPYt28f/vu//xvvv/8+7rvvPs81y5Ytw2uvvYa33noLe/fuxS9+8Qs4HA7cfffdg319REREFGZ8njbatWsXbrzxRs/XSuPsokWLUFZWhsbGRk+QAYCMjAx8/PHHuO+++/Dyyy9j5MiR+POf/4z8/HzPNXfccQdOnjyJFStWwGKx4IorrsDGjRsvauIlIiIiGtQ+L6FioOvEiYiISJyB/v4OnXVPRERERP3A8EJERESqwvBCREREqsLwQkRERKrC8EJERESqwvBCREREqsLwQkRERKrC8EJERESqospTpb9L2WfPZrMJroSIiIj6S/m97et+uWERXlpaWgAA6enpgishIiIiX7W0tMBgMPT7+rA4HsDtduPEiRPQ6/WQJMmvz22z2ZCeno5jx45F5NEDkf76Af4MIv31A/wZRPrrB/gzCNTrl2UZLS0tMJvN0Gj638kSFiMvGo0GI0eODOifkZSUFJF/YRWR/voB/gwi/fUD/BlE+usH+DMIxOv3ZcRFwYZdIiIiUhWGFyIiIlIV7RNPPPGE6CJCnVarxezZsxEVFRazbD6L9NcP8GcQ6a8f4M8g0l8/wJ9BKL3+sGjYJSIiosjBaSMiIiJSFYYXIiIiUhWGFyIiIlIVhhciIiJSFYYXIiIiUhWGl8tYuXIlxowZg9jYWOTm5mLnzp2iSwqK0tJSXH311dDr9UhJSUFBQQH2798vuixhnnvuOUiShOLiYtGlBFVDQwN+/OMfY9iwYYiLi8OUKVOwa9cu0WUFhcvlwvLly5GRkYG4uDiMGzcOTz31lM8HyKnJF198gfnz58NsNkOSJKxfv97r+7IsY8WKFTCZTIiLi0NeXh4OHjwoqFr/u9Tr7+zsxMMPP4wpU6YgISEBZrMZhYWFOHHihMCK/e9yfwcu9POf/xySJOGll14KYoXdGF4u4b333sOyZcvw+OOPY8+ePZg2bRry8/PR3NwsurSAq6ysRFFREbZv347Nmzejs7MTc+bMgcPhEF1a0P373//Gn/70J0ydOlV0KUF19uxZzJo1C9HR0fj0009RW1uLF198EUOHDhVdWlA8//zzWLVqFf74xz9i7969eP755/HCCy/gD3/4g+jSAsbhcGDatGlYuXJlr99/4YUX8Morr2D16tXYsWMHEhISkJ+fj/b29iBXGhiXev2tra3Ys2cPli9fjj179uCvf/0r9u/fjx/84AcCKg2cy/0dUKxbtw7bt2+H2WwOUmXfIVOfrrnmGrmoqMjztcvlks1ms1xaWiqwKjGam5tlAHJlZaXoUoKqpaVFnjBhgrx582b5+9//vnzvvfeKLiloHn74Yfm6664TXYYw8+bNk3/yk5943fef//mf8p133imoouACIK9bt87ztdvtlo1Go/yb3/zGc9+5c+dknU4nv/vuuyJKDKjvvv7e7Ny5UwYgHz16NEhVBVdfP4Pjx4/LaWlpcnV1tTx69Gj597//fdBr48hLHzo6OrB7927k5eV57tNoNMjLy8O2bdsEViaG1WoFACQnJwuuJLiKioowb948r78HkeJvf/sbpk+fjttuuw0pKSm48sor8dprr4kuK2hmzpyJ8vJyHDhwAADwf//3f9i6dSvmzp0ruDIx6urqYLFYvP4tGAwG5ObmRuR7ItD9vihJEoYMGSK6lKBxu92466678OCDDyI7O1tYHeL3+A1Rp06dgsvlQmpqqtf9qamp2Ldvn6CqxHC73SguLsasWbOQk5MjupygWbt2Lfbs2YN///vfoksR4vDhw1i1ahWWLVuGRx99FP/+97/xq1/9CjExMVi0aJHo8gLukUcegc1mQ2ZmJrRaLVwuF5555hnceeedoksTwmKxAECv74nK9yJJe3s7Hn74YSxcuDCiTpl+/vnnERUVhV/96ldC62B4ocsqKipCdXU1tm7dKrqUoDl27BjuvfdebN68GbGxsaLLEcLtdmP69Ol49tlnAQBXXnklqqursXr16ogIL++//z7+8pe/YM2aNcjOzkZVVRWKi4thNpsj4vVT3zo7O3H77bdDlmWsWrVKdDlBs3v3brz88svYs2cPJEkSWgunjfowfPhwaLVaNDU1ed3f1NQEo9EoqKrgW7p0KTZs2IAtW7Zg5MiRossJmt27d6O5uRnf+973EBUVhaioKFRWVuKVV15BVFQUXC6X6BIDzmQyISsry+u+yZMno76+XlBFwfXggw/ikUcewYIFCzBlyhTcdddduO+++1BaWiq6NCGU971If09UgsvRo0exefPmiBp1+ec//zMyBjUAAALCSURBVInm5maMGjXK87549OhR3H///RgzZkxQa2F46UNMTAyuuuoqlJeXe+5zu90oLy/HjBkzBFYWHLIsY+nSpVi3bh0+//xzZGRkiC4pqG666SZ88803qKqq8tymT5+OO++8E1VVVdBqtaJLDLhZs2ZdtDz+wIEDGD16tKCKgqu1tRUajfdbpFarhdvtFlSRWBkZGTAajV7viTabDTt27IiI90SgJ7gcPHgQ//jHPzBs2DDRJQXVXXfdha+//trrfdFsNuPBBx/Epk2bgloLp40uYdmyZVi0aBGmT5+Oa665Bi+99BIcDgfuvvtu0aUFXFFREdasWYOPPvoIer3eM6dtMBgQFxcnuLrA0+v1F/X3JCQkYNiwYRHT93Pfffdh5syZePbZZ3H77bdj586dePXVV/Hqq6+KLi0o5s+fj2eeeQajRo1CdnY2vvrqK/zud7/DT37yE9GlBYzdbsehQ4c8X9fV1aGqqgrJyckYNWoUiouL8fTTT2PChAnIyMjA8uXLYTabUVBQILBq/7nU6zeZTPjRj36EPXv2YMOGDXC5XJ73xeTkZMTExIgq268u93fgu4EtOjoaRqMRkyZNCm6hQV/fpDJ/+MMf5FGjRskxMTHyNddcI2/fvl10SUEBoNfbm2++Kbo0YSJtqbQsy/Lf//53OScnR9bpdHJmZqb86quvii4paGw2m3zvvffKo0aNkmNjY+WxY8fK//Vf/yU7nU7RpQXMli1bev13v2jRIlmWu5dLL1++XE5NTZV1Op180003yfv37xdbtB9d6vXX1dX1+b64ZcsW0aX7zeX+DnyXqKXSkiyH8XaRREREFHbY80JERESqwvBCREREqsLwQkRERKrC8EJERESqwvBCREREqsLwQkRERKrC8EJERESqwvBCREREqsLwQkRERKrC8EJERESqwvBCREREqvL/A97D37+upMYQAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f9d53d488d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x7f9d37776490>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BayesHmm (generic function with 2 methods)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@model BayesHmm(y) = begin\n",
    "    s = tzeros(Int, N)\n",
    "    m = Vector{Real}(K)\n",
    "    T = Vector{Vector{Real}}(K)\n",
    "    for i = 1:K\n",
    "        T[i] ~ Dirichlet(ones(K)/K)\n",
    "        # m[i] ~ Normal(1, 0.1) # Defining m this way causes label-switching problem.\n",
    "        m[i] ~ Normal(i, 0.01)\n",
    "    end\n",
    "    s[1] ~ Categorical(ones(Float64, K)/K)\n",
    "    for i = 2:N\n",
    "        s[i] ~ Categorical(vec(T[s[i-1]]))\n",
    "        y[i] ~ Normal(m[s[i]], 0.01)\n",
    "    end\n",
    "    return(s, m)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Turing]:  Assume - `T` is a parameter\n",
      " @~(::ANY, ::ANY) at compiler.jl:76\n",
      "[Turing]:  Assume - `m` is a parameter (ignoring `m` found in global scope)\n",
      " @~(::ANY, ::ANY) at compiler.jl:76\n",
      "[Turing]:  Assume - `s` is a parameter (ignoring `s` found in global scope)\n",
      " @~(::ANY, ::ANY) at compiler.jl:76\n",
      "[Turing]:  Observe - `y` is an observation\n",
      " @~(::ANY, ::ANY) at compiler.jl:57\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m[Gibbs] Sampling... 99%  ETA: 0:00:01\u001b[39m"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Gibbs] Finished with\n",
      "  Running time    = 53.564168918000036;\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "\u001b[32m[Gibbs] Sampling...100% Time: 0:00:54\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "g = Gibbs(500, HMC(1, 0.2, 3, :m, :T), PG(50, 1, :s))\n",
    "c = sample(BayesHmm(y), g);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iterations = 1:500\n",
      "Thinning interval = 1\n",
      "Chains = 1\n",
      "Samples per chain = 500\n",
      "\n",
      "Empirical Posterior Estimates:\n",
      "               Mean                     SD                           Naive SE                        MCSE                 ESS   \n",
      "T[2][1]  4.792786378×10⁻¹    0.0000000000000015003021363   0.000000000000000067095551270   0.000000000000000000000000 500.000000\n",
      "T[2][2]  2.784222581×10⁻¹    0.0000000000000001667002374   0.000000000000000007455061252   0.000000000000000000000000 500.000000\n",
      "T[2][3]  2.422991041×10⁻¹    0.0000000000000003611838476   0.000000000000000016152632713   0.000000000000000013877788 500.000000\n",
      " lf_num     2.9940000×10⁰    0.1341640786499871118575555   0.005999999999999987981835758   0.005999999999999961093622 500.000000\n",
      "   s[2]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   s[4]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   s[9]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   s[1]     2.0400000×10⁰    0.5127629585914010856839695   0.022931456635085653572581066   0.121202310209005501007162  17.898288\n",
      "   s[6]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "  s[14]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "T[1][1]  8.065717592×10⁻¹    0.0000000000000011113349158   0.000000000000000049700408348   0.000000000000000000000000 500.000000\n",
      "T[1][2]  1.795871242×10⁻¹    0.0000000000000001389168645   0.000000000000000006212551044   0.000000000000000000000000 500.000000\n",
      "T[1][3] 1.3841116588×10⁻²    0.0000000000000000017364608   0.000000000000000000077656888   0.000000000000000000000000 500.000000\n",
      "elapsed  1.071283378×10⁻¹    0.0597890418667208489722498   0.002673847238471375054730261   0.003990129674620669468499 224.527307\n",
      "  s[12]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   s[8]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "  s[11]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "T[3][1]  6.169442702×10⁻¹    0.0000000000000004445339663   0.000000000000000019880163339   0.000000000000000000000000 500.000000\n",
      "T[3][2]  3.444090352×10⁻¹    0.0000000000000004445339663   0.000000000000000019880163339   0.000000000000000000000000 500.000000\n",
      "T[3][3] 3.8646694554×10⁻²    0.0000000000000000138916864   0.000000000000000000621255104   0.000000000000000000000000 500.000000\n",
      "epsilon    2.0000000×10⁻¹    0.0000000000000004723173392   0.000000000000000021122673548   0.000000000000000013877788 500.000000\n",
      "   s[7]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "  s[10]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   m[1]    2.71155628×10⁰    0.0000000000000008890679326   0.000000000000000039760326679   0.000000000000000000000000 500.000000\n",
      "   s[3]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "   m[2]    1.89621886×10⁰    0.0000000000000015558688821   0.000000000000000069580571688   0.000000000000000111022302 196.392786\n",
      "   s[5]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "  s[15]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "  s[13]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN\n",
      "     lp   -4.46911188×10⁴ 3656.1798791090955091931391507 163.509334953098033338392269798 163.360946325820577840204351 500.000000\n",
      "   m[3]  7.111201435×10⁻¹    0.0000000000000010002014242   0.000000000000000044730367514   0.000000000000000000000000 500.000000\n",
      "\n",
      "Quantiles:\n",
      "               2.5%             25.0%             50.0%             75.0%             97.5%      \n",
      "T[2][1]  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹\n",
      "T[2][2]  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹\n",
      "T[2][3]  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹\n",
      " lf_num     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "   s[2]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "   s[4]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "   s[9]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰\n",
      "   s[1]     1.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     3.0000000×10⁰\n",
      "   s[6]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "  s[14]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "T[1][1]  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹\n",
      "T[1][2]  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹\n",
      "T[1][3] 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻²\n",
      "elapsed 8.2674258675×10⁻²   9.60775315×10⁻²   1.06401758×10⁻¹  1.098029557×10⁻¹  1.302064362×10⁻¹\n",
      "  s[12]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "   s[8]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰\n",
      "  s[11]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "T[3][1]  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹\n",
      "T[3][2]  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹\n",
      "T[3][3] 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻²\n",
      "epsilon    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹\n",
      "   s[7]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "  s[10]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰\n",
      "   m[1]    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰\n",
      "   s[3]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "   m[2]    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰\n",
      "   s[5]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "  s[15]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰\n",
      "  s[13]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰\n",
      "     lp   -4.45298711×10⁴   -4.45288443×10⁴   -4.45270086×10⁴   -4.45270086×10⁴   -4.45270086×10⁴\n",
      "   m[3]  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹\n",
      "\n"
     ]
    }
   ],
   "source": [
    "describe(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGgCAYAAABxDccgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xt0VPW9///XJCGXahIEzU0SCBbDHbmIAlqooJSvXy2r69hf/VnB2vptz4FvRbo81p6j2ONRtGf1eKkWql2VquXQ2hbt8VRtjIBS5RoCCUK4iFwTLqJJiJDEzP7+8XHnIknIJDP7s2fm+Vhr1uwMM3zemezZ+zWf/dmfHXAcxxEAAIBPJNguAAAAoC3CCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8JUk2wV0RzAY1JEjR5Senq5AIGC7HAAA0A2O46iurk55eXlKSOh+f0hUhJMjR44oPz/fdhkAAKAHDh48qAEDBnT7+VERTtLT0yWZXy4jI8NyNQAAoDtqa2uVn5/fsh/vrqgIJ+6hnIyMDMIJAABRJtQhGQyIBQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AQAAvkI4AWBXVZX0wAPmnvYBKMRwsmTJEo0ePbplGvlJkybptdde6/I1L730koYOHarU1FSNGjVKf/3rX3tVMIAYU1Ul/fSndsNBPLcP+FBI4WTAgAF65JFHtHnzZm3atEnXXHONvv71r2v79u0dPv/dd9/VzTffrO9+97vasmWLZs+erdmzZ6uioiIsxQMAgNgT0oX/brjhhnY/P/TQQ1qyZInWrVunESNGnPX8J554Ql/72td09913S5IefPBBFRcX66mnntLSpUs7baehoUENDQ0tP9fW1oZSJgC/q6pq7SkoLW1/L0m5ueZG+0Bc6vGYk+bmZq1YsUL19fWaNGlSh8957733NGPGjHaPzZw5U++9916X//fixYuVmZnZcsvPz+9pmQD86Fe/ksaPN7c77jCP3XFH62O/+hXtA3EspJ4TSSovL9ekSZN05swZnX/++Vq5cqWGDx/e4XOrq6uVnZ3d7rHs7GxVV1d32ca9996rhQsXtvxcW1tLQAFiyfe/L914o1kuLTU75meflcaNM49Futcg3tsHfC7kcFJUVKSysjLV1NToj3/8o+bOnas1a9Z0GlB6IiUlRSkpKWH7/wD4TEeHLcaNa9050z4Q10IOJ8nJyfryl78sSRo/frw2btyoJ554Qr/qoBsyJydHR48ebffY0aNHlZOT08NyAQBArOv1PCfBYLDd4NW2Jk2apJKSknaPFRcXdzpGBUAcys2VFi2ydygj3tsHfCjgOI7T3Sffe++9mjVrlgoKClRXV6fly5fr0Ucf1RtvvKFrr71Wc+bM0cUXX6zFixdLMqcST506VY888oiuv/56rVixQg8//LBKS0s1cuTIbhdZW1urzMxM1dTUKCMjI/TfEgAAeK6n+++QDuscO3ZMc+bMUVVVlTIzMzV69OiWYCJJBw4cUEJCa2fM5MmTtXz5cv3rv/6rfvKTn2jIkCF6+eWXQwomAAAgvoTUc2ILPScAAESfnu6/ubYOAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwFcIJAKs2bZKuucbc0z4AiXACwLLnn5dWrZJeeIH2ARhJtgsAEH/275dOnJACAWn5cvPYCy9IY8dKjiNlZkq5uZFrv6pKqqkx7b/4oj/aX7FCmjvXtH/hhdLAgZFrH/C7gOM4ju0izqW2tlaZmZmqqalRRkaG7XIA9FIgYLsC/wkETDBx+X/LDJxbT/ff9JwA8NyLL0q33SZ99lnH/37RRVJ6euTar6uTjh/v/N9ttO+GkaQkadmyyLUNRAPCCQDP3XKLNGyYNH782f+2ebM0blzkaygt9Wf769d70z7gZwyIBWDFRx+1/znB0tbIbddW+xziAs5GOAFgxc6d5j4tTVq61PQi5ORIWVnetJ+VZdobP95++6mp5rH+/b1rH/AzDusAsKK01NzPmyd9//vS//k/UmOjlJLiTfsDBkgffiglJ5veC5vtf+Mb0ssvSz/8oXkciHf0nADwnONIxcVm+dprzX0g4F0wcKWktB5Wsdn+jBnm51WrvG0f8CvCCQDP7dwpHT5sds5XX227GvvcgPb3v0v19XZrAfyAcALAc26vyVVXmTEn8W7IEKmgQGpqkt5+23Y1gH2EEwCee/NNc+/2GMS7QKD1vXDfGyCeEU4AeKqpSVq92iwTTlq574XbqwTEM8IJAE+tX29mSL3wQumyy2xX4x/Tp5selPJyqbradjWAXYQTAJ5yewamT7c38ZkfXXihufCgxKEdgE0DAE998RRitOLQDmAQTgB4pqZG2rDBLBNOztY2nHBVYsQzwgkAz6xaJTU3t546i/amTDFT2VdVSe+/b7sawB7CCQDPcEina6mprZPScWgH8YxwAsAzhJNzY9wJQDgB4JH9+6Xdu6XEROmrX7VdjX+54WTNGnMhQiAeEU4AeMI9PXbiRCkz024tfjZ6tHTRReYaO+vW2a4GsINwAsATHNLpnoSE1qsUc2gH8YpwAiDigkGppMQsE07OjXEniHeEEwARV1YmnTghpadLV1xhuxr/c8PJxo3Sxx/brQWwgXACIOLcHoBp06Q+fayWEhUGDJCGDjU9TqtW2a4G8B7hBEDEueHEHUuBc2PcCeIZ4QRARJ0+La1da5YZb9J9jDtBPCOcAIiod96RGhqkiy82hyrQPdOmmTlh9u6V9u2zXQ3gLcIJgIhy5ze59lopELBbSzTJyJCuvNIsu+8hEC9CCieLFy/W5ZdfrvT0dGVlZWn27NmqrKzs8jXLli1TIBBod0tNTe1V0QCiB/Ob9ByHdhCvQgona9as0bx587Ru3ToVFxerqalJ1113nerr67t8XUZGhqqqqlpu+/fv71XRAKLDsWPmNGKJwbA94YaTkhJzNWcgXiSF8uTXX3+93c/Lli1TVlaWNm/erK985Sudvi4QCCgnJ6dnFQKIWu7Ea2PGSFlZdmuJRhMnmsM7J09KW7ZIEybYrgjwRq/GnNTU1EiS+vXr1+XzTp06pYEDByo/P19f//rXtX379i6f39DQoNra2nY3ANGHQzq9k5TUepFEDu0gnvQ4nASDQS1YsEBTpkzRyJEjO31eUVGRfvOb3+iVV17Riy++qGAwqMmTJ+vQoUOdvmbx4sXKzMxsueXn5/e0TACWOA7hJBwYd4J4FHAcx+nJC//xH/9Rr732mtauXasBAwZ0+3VNTU0aNmyYbr75Zj344IMdPqehoUENDQ0tP9fW1io/P181NTXKyMjoSbkAPLZzpzRsmJScbKZg/9KXbFcUnSorzSnYvI+IRrW1tcrMzAx5/92jnpP58+fr1Vdf1apVq0IKJpLUp08fjR07Vnv27On0OSkpKcrIyGh3AxBd3G/6V13FDrU3Lr1Uys+XGhult9+2XQ3gjZDCieM4mj9/vlauXKm33npLhYWFITfY3Nys8vJy5ebmhvxaANGDQzrhEQhwaAfxJ6RwMm/ePL344otavny50tPTVV1drerqap0+fbrlOXPmzNG9997b8vO//du/6W9/+5s++OADlZaW6tvf/rb279+v733ve+H7LQD4SlOTtHq1WSac9J77HjIZG+JFSKcSL1myRJI0bdq0do8/99xzuu222yRJBw4cUEJCa+b5+OOPdccdd6i6uloXXHCBxo8fr3fffVfDhw/vXeUAfGvDBqmuTurfXxo71nY10W/6dHO/bZt09KiUnW23HiDSQgon3Rk7u9r9uvS5xx57TI899lhIRQGIbu7hh+nTpQQuktFrF11kQt6WLab35JZbbFcERBabDQBhx3iT8GPcCeIJ4QRAWNXUSOvXm2XCSfi0DSc9mwACiB6EEwBhtXq1uQ7MkCHSwIG2q4kdV10lpaZKR45IO3bYrgaILMIJgLByDztwob/wSk01AUXi0A5iH+EEQFgx3iRyGHeCeEE4ARA2Bw5Iu3aZM3TcC9YhfNxwsnq1mTEWiFWEEwBh404SNnGi1Lev3Vpi0Zgx5rTi+vrWQcdALCKcAAgbDulEVkJC64RsHNpBLCOcAAiLYLC154RwEjmMO0E8IJwACIutW6UTJ6Tzz5euvNJ2NbHLDScbNkiffGK3FiBSCCcAwsL9Jj9tmtSnj9VSYlp+vlRUZHqqVq2yXQ0QGYQTAGHBeBPvcGgHsY5wAqDXTp+W3nnHLDP5WuS57zHhBLGKcAKg19aulRoapLw8adgw29XEvmnTpMREac8e6cMPbVcDhB/hBECvtT2kEwjYrSUeZGZKV1xhluk9QSwinADoNU4h9p77XrvvPRBLCCcAeuX4cWnLFrPMeBPvuOGkpMScuQPEEsIJgF4pKTH3o0dL2dl2a4knEydK6enSRx+1hkMgVhBOAPQKpxDb0adP68UVGXeCWEM4AdBjjkM4sYn5ThCrCCcAemzXLungQSk5Wbr6atvVxB83nKxdK336qd1agHAinADoMfcb+5Qp0pe+ZLeWeHTppdKAAVJjY+skeEAsIJwA6DEO6dgVCHBoB7GJcAKgRz77rPXCc4QTe5jvBLGIcAKgRzZskOrqpH79pLFjbVcTv6ZPN/dbt0pHj9qtBQgXwgmAHnEPI0yfbq7zAjuysqTLLjPL7pwzQLQjnADoEcab+AfjThBrCCcAQlZbK61bZ5YJJ/a1DSeOY7cWIBwIJwBCtnq11NwsffnL0qBBtqvBVVdJKSnS4cPSzp22qwF6j3ACIGQc0vGXtLTWSfA4tINYQDgBEDJ3B8hViP3D/VsQThALCCcAQnLwoFRZKSUkSNdcY7sauNxerNWrpaYmq6UAvUY4ARASd7Kvyy+X+va1WwtaXXaZdOGF0qlT0vr1tqsBeodwAiAkjDfxp4SE1gnZOLSDaEc4AdBtwWBrzwnhxH+Y7wSxgnACoNu2bZOOH5fOO0+68krb1eCL3HCyYYNUU2O3FqA3CCcAus39Rj5tmpScbLUUdKCgQLr0UjMHjXtRRiAaEU4AdBvjTfyPQzuIBYQTAN1y5oz0zjtmmXDiX4QTxALCCYBuWbvWBJS8PGnYMNvVoDPTppmrRO/eLe3fb7saoGcIJwC6pe2ssIGA3VrQucxMaeJEs0zvCaIV4QRAt3AKcfRw/0bu3wyINoQTAOd04oS0ZYtZ5no6/ueGk5ISMzcNEG0IJwDOqaREchxp1CgpJ8d2NTiXK66Q0tNNqCwrs10NEDrCCYBz4hTi6NKnjxkYKzHuBNGJcAKgS45DOIlGnFKMaEY4AdCl3bulAwfMjLBXX227GnSXG07WrpVOn7ZbCxAqwgmALrnfvCdPNtfUQXQoKpIuvlhqaGidPA+IFoQTAF3ikE50CgQ4tIPoRTgB0KnPPmu9gBzhJPoQThCtCCcAOrVxo1RbK11wgTRunO1qECp3TpqtW6Vjx+zWAoSCcAKgU+437unTzfVaEF2ysqQxY8xySYndWoBQEE4AdIrxJtGPQzuIRoQTAB2qq5PWrTPLhJPo1TacOI7dWoDuIpwA6NDq1WZA7CWXSIWFtqtBT119tZSSIh06JFVW2q4G6B7CCYAOcUgnNqSlSVddZZY5tINoQTgB0CF3R8ZViKOf+zcknCBaEE4AnOXQIWnnTikhQbrmGtvVoLfc3q/Vq6WmJqulAN1COAFwFvcb9oQJZo4TRLexY6X+/c0g5/XrbVcDnBvhBMBZ3nzT3DPeJDYkJJi5aqTWvy3gZ4QTAO0Eg4STWMR8J4gmhBMA7ZSXm6nOzztPmjTJdjUIFzecrF8v1dTYrQU4F8IJ4ltVlfTAA+Y+HtvvoAb3m/XUqVJysr2yEF4DB0pDhkjNzWZgbAsfroNx1z7OElI4Wbx4sS6//HKlp6crKytLs2fPVmU3ZvV56aWXNHToUKWmpmrUqFH661//2uOCgbCqqpJ++lO7G0Wb7XdQA/ObxK4OD+34cB2Mu/ZxlpDCyZo1azRv3jytW7dOxcXFampq0nXXXaf6+vpOX/Puu+/q5ptv1ne/+11t2bJFs2fP1uzZs1VRUdHr4gGE15kz0ttvm2XCSexh3AmiRVIoT3799dfb/bxs2TJlZWVp8+bN+spXvtLha5544gl97Wtf09133y1JevDBB1VcXKynnnpKS5cu7fA1DQ0NamhoaPm5trY2lDKBrlVVtX5DKi1tfy9JubnmFqvtd1HD3zek68yZIcrNbtbw4VyGONZMm2bO3Nm1SzrwP+UqyG3y3TroWQ2220fXnF7YvXu3I8kpLy/v9Dn5+fnOY4891u6x+++/3xk9enSnr1m0aJEj6axbTU1Nb8oFjEWLHMdcA63j26JFsd1+FzXco8WO5Di3ji6LfA2w4sorzZ/717rdl+tgXH0O40BNTU2P9t8Bx+nZdSqDwaBuvPFGffLJJ1q7dm2nz0tOTtZvf/tb3XzzzS2P/fKXv9RPf/pTHT16tMPXdNRzkp+fr5qaGmVkZPSkXKDVF78x3XGH9Oyz0rhx5jGvv7F53X4XNYz/+f+v0p1f0vNPfqJb/2/fyNYAK+6/X3rwQen/u+6kViz+0HfrYFx9DuNAbW2tMjMzQ95/h3RYp6158+apoqKiy2DSUykpKUpJSQn7/wtI6nijM25c60Yp1tvvpIYThZdrS+WXJEkz/oFgEquuvdaEk5LSfgpe1q914KEP1sG4+xyiUz06lXj+/Pl69dVXtWrVKg0YMKDL5+bk5JzVQ3L06FHl5OT0pGkAEfLWxnQ5jjRyJF8YY9mVV0rnny+dOCFt3Wq7GqBjIYUTx3E0f/58rVy5Um+99ZYKCwvP+ZpJkyappKSk3WPFxcWaxOxO8IPcXGnRInt7Y9vtt6mhuNx8YeAsndjWp48ZGCt9ftaOj9bBuP4cop2Qxpz80z/9k5YvX65XXnlFRUVFLY9nZmYqLS1NkjRnzhxdfPHFWrx4sSRzKvHUqVP1yCOP6Prrr9eKFSv08MMPq7S0VCNHjuxWuz09ZgWgexxHKiyU9u+X/vpXadYs2xUhkp58UrrzTmnGDE4rRmT1dP8dUs/JkiVLVFNTo2nTpik3N7fl9vvf/77lOQcOHFBVm4lsJk+erOXLl+uZZ57RmDFj9Mc//lEvv/xyt4MJgMjbs8cEk+RkqZNZARBD3N6xd96RTp+2WwvQkZAGxHank2V1u3mRjZtuukk33XRTKE0B8JD77XnyZHNNHcS2oUOliy+WDh+W1q7lUB78h2vrAGgJJzNm2K0D3ggEWv/WHNaBHxFOgDj32WfSW2+ZZb5Bxw+msoefEU6AOLdxo1RbK11wgTR+vO1q4BW356SsTDp2zG4twBcRToA49+ab5v6aa6RELqcTN7KzpdGjzbLbcwb4RVyHk02bzAZ50ybblcAW2+uA7fYl6c9/NveXXmqvBtjhHtq56y6766Dtz4Ht9nG2uA4nzz8vrVolvfCC7Upgi+11wHb7dXWts4QeOWKnBtjjhpPqarMu2mL7c2C7fZytx9fWiVb795tpmwMBadky89jvfifNnWsmorrwQmngQKslIsLargPuxujXv5ZqasxyWpqUmRm59mtqWueWeOkl79tvW8OxY2a9l6TXXjPXP+NzEPvcz8D557c+9utfm7FHkrfroGTnc9BR+ytWsC/wix5fldhL4ZwhNhA493P8/46gN7qzDsSjQKD9us/nIHbxGegYn4Hw8/yqxNHqxRel224zp09+UVJSa28KYldX60BCgjRzpjRsWOTa37FDeuMNKRi0035nNbgbYj4Hsc/2Z0Cy/zngM+BvcddzIpmu645Omdy8matlxwvb64Dt9v1SA+zxw9/fdg22248HnlxbJ9bQtQlXgqVPgtuurfb9UgPs8cPf3w81SOwT/CQuN0dZWVJOjnTZZe0fy8qyVxO8lZUlpaSY5W99y3x7ysnxbh1w18Hx46WlS71v3y81wB4//P1t1+C2n53d/mc+A/bF5WEdSWpoMFdgHTxY+vBDc+zxuuvC8l8jSgwYYC589s470pQpUmNja2DxgrsOuoPwvG7fLzXAHj/8/W3X0NBgztj87neladOk11/nMxBOHNYJUUqK+TC4vSc7dtitB9766CMTTCQzS2Yg4P0GyV0HJTvt+6UG2OOHv7/tGlJSpLFjzfLWrSYowb64DSeuMWPMfVmZ3TrgLXfiscGDpTB1xgGIUsOHmzN0Pv5YOnTIdjWQCCct4cTdWSE+uH/vtuOOAMSnlJTW05bZF/hD3IcTd+e0fbvU1GS3FnjH7SlzwymA+EYvur/EfTgZNMh06zc2Sjt32q4GXnG/HRFOAEj0ovtN3IeTQICVMt40Nkrvv2+WOawDQGrdFtBz4g9xH04kuvPizY4d5hBe375SQYHtagD4gbsf2LtXOnXKbi0gnEii5yTetD2kw4yQACTpooukvDwz10p5ue1qQDhRa3fe1q1chTIeMBgWQEfoRfcPwomkESPMNR2OH5eqqmxXg0hjMCyAjtCL7h+EE0lpadLQoWaZlTK2OU7rtyIGwwJoi0Gx/kE4+RzdefHh8GHp5EkpMdHMCgkALnc/UF4uNTfbrSXeEU4+R3defHD/vsOGSampdmsB4C9Dhpie9E8/NWftwB7CyefaDopF7GIwLIDOJCZKo0aZZXrR7SKcfM7dWe3aZVIzYhODYQF0hV50fyCcfC4nR8rOloJBqaLCdjWIFAbDAugKg2L9gXDSBoNiY1t9vbRnj1mm5wRAR+g58QfCSRuslLGtvNycSpyTI2Vl2a4GgB+NHm3uDx+WTpywW0s8I5y0waDY2MYhHQDnkp4uXXKJWWZfYA/hpI22PSfBoN1aEH4MhgXQHfSi20c4aaOoSEpJMVek3LfPdjUIN3pOAHQHg2LtI5y0kZQkjRxpllkpY0sw2HqlUXpOAHSFnhP7CCdfwEoZm/buNWfrpKaaWSABoDPufmDHDqmx0W4t8Ypw8gUMio1Nbk/YqFGmhwwAOlNQIPXtKzU1Se+/b7ua+EQ4+QLmOolNDIYF0F2BAL3othFOvsBdIQ8ckD7+2G4tCB8GwwIIBYNi7SKcfEFmpjRokFkmMccOek4AhIKeE7sIJx1gpYwtH30kHTpklt3ZHwGgK233A45jt5Z4RDjpAINiY4v7dxw8WMrIsFsLgOgwfLgZPH/yZOuXG3iHcNIBBsXGFg7pAAhVaqo0dKhZ5ouq9wgnHXB7TrZvN6eSIboxGBZATzAo1h7CSQcGDTLd/42N0s6dtqtBb9FzAqAnGH9oD+GkA4FA68BJVsro1tjYOokS4QRAKAgn9hBOOsGg2NiwY4c5NJeZKQ0caLsaANHEDSd79pgLwsI7hJNOMCg2NrQ9pBMI2K0FQHTJypJyc82pxO6FQ+ENwkkn2vaccI579GIwLIDeYFCsHYSTTowYISUkSMePS1VVtqtBTzEYFkBvMO7EDsJJJ9LSpKIis8xKGZ0ch3ACoHcIJ3YQTrrAoNjodviwmbo+MdH0hAFAqNz9wLZtUnOz3VriCeGkCwyKjW5uqBw61Mz2CAChGjLE9KR/+qm0d6/tauIH4aQLdOdFNwbDAuitxERp1CizzBdV7xBOuuDu1HbtMqkZ0YXxJgDCgS+q3iOcdCEnx5znHgxKFRW2q0GoCCcAwoFw4j3CyTkwKDY61ddLu3ebZcIJgN5grhPvEU7OgUGx0am83JxKnJMjZWfbrgZANHOvteaeAYjII5ycA9150cn9ezEYFkBvpadLl1xiltkXeINwcg5tD+sEg3ZrQfe5PV0c0gEQDvSie4twcg5FRVJKirki5b59tqtBdzEYFkA40YvuLcLJOSQlSSNHmmVWyugQDJrZHCUO6wAIDwbFeivkcPL222/rhhtuUF5engKBgF5++eUun7969WoFAoGzbtXV1T0u2mt050WXvXvN2TqpqWZ2RwDoLXc/sGOH1Nhot5Z4EHI4qa+v15gxY/T000+H9LrKykpVVVW13LKyskJt2hq686KL+3caNcr0fAFAbxUUSH37Sk1NJqAgskLedM+aNUuzZs0KuaGsrCz17ds35Nf5Ad150YXBsADCLRAw25Q1a8w2hu1LZHk25uSyyy5Tbm6urr32Wv3973/v8rkNDQ2qra1td7PJPcf9wAHp44+tloJuYDAsgEigF907EQ8nubm5Wrp0qf70pz/pT3/6k/Lz8zVt2jSVlpZ2+prFixcrMzOz5Zafnx/pMrvUt680aJBZdgdawr+44B+ASKAX3TsBx3GcHr84ENDKlSs1e/bskF43depUFRQU6IUXXujw3xsaGtTQ0NDyc21trfLz81VTU6OMjIyeltsrs2dLr7wiPf64dOedVkpAN5w8KfXvb5ZraiRLqwuAGFRaKo0fL/XrJ504YQ71oGu1tbXKzMwMef9t5VTiiRMnas+ePZ3+e0pKijIyMtrdbKM7Lzq4f5/BgwkmAMJr+HAzyP7kSTOVPSLHSjgpKytTbm6ujaZ7jO686MBgWACRkpoqDR1qltkXRFbIZ+ucOnWqXa/Hvn37VFZWpn79+qmgoED33nuvDh8+rOeff16S9Pjjj6uwsFAjRozQmTNn9Otf/1pvvfWW/va3v4Xvt/CAu7Pbvt2cStanj9160DEGwwKIpDFjpIoKs6353//bdjWxK+Sek02bNmns2LEaO3asJGnhwoUaO3as7r//fklSVVWVDhw40PL8xsZG/ehHP9KoUaM0depUbd26VW+++aamT58epl/BG4MGmcMEjY1SZaXtatAZBsMCiCR60b3RqwGxXunpgJpwu/pqae1a6YUXpG9/21oZ6ERjo3T++aZna9++1jOsACBcioul664zs0/5NcnYAAAX3UlEQVTv2mW7Gv+LqgGx0YpBsf62c6cJJpmZ0sCBtqsBEIvc/cCePeaCsIgMwkkI6M7zt7aDYTnFD0AkZGVJubmS40jl5bariV2EkxC07Tnx/8Gw+MNgWABeoBc98ggnIRg5UkpIkI4fl6Loospxg8GwALxAL3rkEU5CkJYmFRWZZVZKf3Ecek4AeIOek8gjnISIldKfjhyRPvpISkyURoywXQ2AWObuB8rLpeZmu7XEKsJJiOjO8yf37zF0qJnFEQAi5dJLTU96fb20d6/tamIT4SRE9Jz4E4d0AHglMdGMQZTYF0QK4SREbs/Jrl3Sp5/arQWtGAwLwEv0okcW4SREOTnmPPdg0FxfAf5AzwkAL9GLHlmEkx5gpfSX+npp926zTDgB4AX2A5FFOOkBuvP8pbzcnEqckyNlZ9uuBkA8GD3a3B86ZM4URHgRTnqAxOwvHNIB4LWMDGnwYLPMviD8CCc94PacbNtmxp7ALgbDArCBXvTIIZz0QFGRlJIi1dVJ+/bZrgb0nACwgV70yCGc9EBSUusspKyUdgWDpgdLIpwA8BbhJHIIJz1Ed54/7N1rztZJTTWzNgKAV9z9wPvvS42NdmuJNYSTHiIx+4P7/o8caXq0AMArBQVS375SU5O0Y4ftamIL4aSH3MRMOLGLwbAAbAkEWr+o0oseXoSTHnLPcd+/X/r4Y7u1xDMGwwKwiV70yCCc9FDfvtLAgWbZHZAJ7xFOANhEOIkMwkkvcGjHrpMnpYMHzbLbkwUAXmp7coTj2K0llhBOeoFjjXa5obCwUMrMtFsLgPg0fLiUmGi+LB0+bLua2EE46QV6TuxiMCwA21JTpWHDzDJfVMOHcNILbs9JRYU5lQzeYrwJAD9g3En4EU56YdAgKT3dTL5TWWm7mvhDOAHgB4ST8COc9EJCAiulLY2N0vbtZpnDOgBsYsbw8COc9BKDYu3YudMcSsvMbD2lGwBscPcDe/ZIp07ZrSVWEE56iUGxdrhhcMwYM0sjANiSlSXl5ppTicvLbVcTGwgnvdS254Rz3L3DeBMAfsIh/vAinPTSyJFm7Mnx41J1te1q4gfhBICfEE7Ci3DSS2lpUlGRWWal9IbjMMcJAH9hUGx4EU7CgEGx3jpyRProIzMr44gRtqsBgNb9QHm5FAzarSUWEE7CgEGx3nLf56FDzeyMAGDbpZeanvT6emnvXtvVRD/CSRjQc+KttmfqAIAfJCaaMYgS+4JwIJyEgbuT3LVLOn3abi3xgMGwAPyIQbHhQzgJg5wcc557MGius4PIYjAsAD9iUGz4EE7CIBDg0I5X6uul3bvNMj0nAPyEnpPwIZyECSulNyoqzKnEOTlSdrbtagCg1ejR5v7QIXNGIXqOcBImdOd5g8GwAPwqI0MaPNgs80W1dwgnYeLuLLdt4xz3SGIwLAA/oxc9PAgnYVJUJKWkSHV10ocf2q4mdjEYFoCf0YseHoSTMOnTp3W2UlbKyAgGTc+URM8JAH+i5yQ8CCdhxEoZWR98YM7WSU01szECgN+4PSfvvy81NtqtJZoRTsKI7rzIct/XkSOlpCS7tQBARwoKpL59paYmaccO29VEL8JJGNFzElkMhgXgd4FA6ynF7At6jnASRu5Oc/9+6ZNP7NYSixgMCyAa0Ivee4STMOrbVxo40CyTmMOPnhMA0YBe9N4jnIQZK2VknDwpHTxolt0uUwDwo7b7AcexW0u0IpyEGd15keGGvcJCKTPTbi0A0JURI6TERDOF/eHDtquJToSTMKPnJDI4pAMgWqSmSkOHmmX2BT1DOAkzt+dk+3ZzKhnCg8GwAKIJvei9QzgJs0GDpPR0qaFBqqy0XU3soOcEQDShF713CCdhlpDAOe7h1thoZluUCCcAogPhpHcIJxFAd1547dxpAkpGhumZAgC/c8PJ7t3mshsIDeEkAkjM4dX2kE4gYLcWAOiO7GwpJ8ecSlxebrua6EM4iYC2PSec4957DIYFEI3oRe85wkkEjBxpxp4cPy5VV9uuJvoxGBZANKIXvecIJxGQliZdeqlZZqXsHcchnACIToSTniOcRAjdeeFx5Ih04oSZbXHECNvVAED3ufuBbdukYNBuLdGGcBIhJObwcN+/oiLTIwUA0WLIEDNbbH29tHev7WqiC+EkQtzETDjpHQbDAohWSUnSqFFmmV700IQcTt5++23dcMMNysvLUyAQ0Msvv3zO16xevVrjxo1TSkqKvvzlL2vZsmU9qTX8qqqkBx4w92Hm9pxUVkqnT3vffrfZruEc7TPeBEA061Yvus+3wzaEHE7q6+s1ZswYPf300916/r59+3T99dfrq1/9qsrKyrRgwQJ973vf0xtvvBFysWFXVSX99KcR+YPk5EgXXWSOM1ZUeN9+t9mu4RztE04ARLNuhxMfb4dtSAr1BbNmzdKsWbO6/fylS5eqsLBQP//5zyVJw4YN09q1a/XYY49p5syZoTYfNQIBcyiiuNh0511+ue2Kok99vbRrl1nmsA6AaMTJET0TcjgJ1XvvvacZM2a0e2zmzJlasGBBp69paGhQQ0NDy8+1tbXhK6iqqjUdlpa2v5ek3FxzC4MxY0w4aZeYPWy/U7Zr6Gb7FRXmVOLsbHMDgGjjXmvt0CHpo4+k/v0//4co2Q5b4/SCJGflypVdPmfIkCHOww8/3O6x//mf/3EkOZ9++mmHr1m0aJEj6axbTU1Nb8p1/3PHMfu8jm+LFvW+jc+9+KL5L6+6yk77nbJdQzfbX7rU/DhzZmTLAYBIGjzYbMtKSto8GCXb4d6qqanp0f474j0nPXHvvfdq4cKFLT/X1tYqPz8/PP/5978v3XijWS4tle64Q3r2WWncOPNYGJNi22ONwaCZNdbL9jtlu4Zuts94EwCxYMwY6YMPzDbtmms+fzBKtsO2RDyc5OTk6OjRo+0eO3r0qDIyMpTWycQVKSkpSklJiUxBHXVVjRvX+gcJo6IiKTlZqquTPvxQGjzY2/Y7ZbuGbrZPOAEQC8aMkVau/MIh/ijZDtsS8XlOJk2apJKSknaPFRcXa9KkSZFu2ro+fcx1diQGQ4UqGGz9IDMYFkA0Y1Bs6EIOJ6dOnVJZWZnKPn+X9+3bp7KyMh04cECSOSQzZ86cluf/4Ac/0AcffKB//ud/1s6dO/XLX/5Sf/jDH3TXXXeF6VfohdxcadGiiHZfdXkamQftn5PtGjpp/4MPzNk6KSmt1ykCgGjk7gfef19qbOzgCT7dDlsV6uCWVatWdThYde7cuY7jOM7cuXOdqVOnnvWayy67zElOTnYGDx7sPPfccyG12dMBNX7wxBNmbNHXv267kujy0kvmfZswwXYlANA7waDjZGaabVpZme1qvOXZgNhp06bJcZxO/72j2V+nTZumLVu2hNpUTHATM915oWG8CYBYEQiYbdnbb5ttG9u1c+PaOhHmnuO+f7/0ySd2a4kmhBMAsYSLwYaGcBJhF1wgDRxollkpu48L/gGIJQyKDQ3hxAMk5tCcPCkdPGiW3Z4nAIhmbfcDXYyMwOcIJx5wEzPhpHu2bTP3hYVSZqbdWgAgHEaMkBITzRT2R47Yrsb/CCceYFBsaNz3ifEmAGJFaqo0dKhZZl9wboQTD7g72e3bpaYmu7VEAwbDAohFHOLvPsKJBwoLpfR0qaFBqqy0XY3/MRgWQCxiUGz3EU48kJDQOrCTxNy1piYzi6JEzwmA2ELPSfcRTjzCoNju2bnTTO+ckSENGmS7GgAIHzec7N5tLs+BzhFOPMKg2O5pOxg2ELBbCwCEU3a2lJNjTiUuL7ddjb8RTjzSNpxwjnvnGAwLIJZxaKd7CCceGTnSjD05flyqrrZdjX8xGBZALGNQbPcQTjzypS9Jl15qlknMHXMcek4AxDZ6TrqHcOIhBsV2rapKOnHCzKI4YoTtagAg/Nz9wLZtUjBotxY/I5x4iEGxXXPfl6IiKS3Nbi0AEAlDhpjZYuvrpb17bVfjX4QTD9Gd1zUO6QCIdUlJZgyixL6gK4QTD7ndeZWV0unTdmvxIwbDAogHDIo9N8KJh3JypIsuMscZKypsV+M/9JwAiAf0op8b4cRDgQArZWfq66Vdu8wyPScAYhknR5wb4cRjdOd1rKLCnEqcnW1uABCr3GutHTwonTxptxa/Ipx4jJ6TjnFIB0C8yMgwV6uX2Bd0hnDisbbdeZzj3orBsADiCb3oXSOceKyoSEpOlurqpA8/tF2Nf9BzAiCe0IveNcKJx/r0aZ39lJXSCAbNbIkS4QRAfCCcdI1wYgHdee198IF06pSUkmJ6lgAg1rn7ge3bpcZGu7X4EeHEAhJze+77MHKkmT0RAGLdwIFSZqbU1CTt3Gm7Gv8hnFhAz0l7DIYFEG/aznvFvuBshBML3HPc9++XPvnEbi1+wGBYAPGIXvTOEU4suOACqaDALLsDQeMZ4QRAPCKcdI5wYgmHdoyTJ6UDB8wy4QRAPGm7H3Acu7X4DeHEEhKz4fYcDRpkBocBQLwYMUJKTJQ++kg6csR2Nf5COLHETcx/+IO0aZO9OjZtkq65xl4Nf/mLuR840E77AGBLaqo0dKhZ/l//y9522PZ+oCOEE0vcnpNTp6Tf/tZeHc8/L61aJb3wgp32X3/d3J8+bad9ALDJ3Rds22ZvO2x7P9ARZpXw2P790okT7Y8vLl8ufec75rELL4x8L4JbQyAg/f735rEVK6S5c72poW37u3ebxyorpdJS794DALDJ3Q5eeGHrY7a2wzb2A+cScBz/D8Opra1VZmamampqlJGRYbucXgkEzv2cO++MbA1PPGG3hu607/+1EgB6zva+oKPtcCDQftsbju1wT/ff9Jx47MUXpdtukz77rPPndGfnHWm2akhKkpYts9M2AHjFj/sCN4z4YTtMOPHYLbdIw4ZJ48ef/W+33y7l5HhTR3W19Jvf2Kuhs/bXr5fGjYt8+wBgkx/2BX7eDhNOLEpIMFfkde/nzfNuhSgtNSulrRo6ax8A4g3b4bNxto4FWVkmFY8fLy1dau5zcszj8VKD7fYBwDbb20Hb7XeFAbGWNDRIycmtA5AaG6WUlPiqwXb7AGCb7e1gpNtnQGyUafvHDwTs7JRt12C7fQCwzfZ20Hb7neGwDgAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8BXCCQAA8JWomL7evfxPbW2t5UoAAEB3ufvtUC/jFxXhpK6uTpKUn59vuRIAABCquro6ZWZmdvv5UXFV4mAwqCNHjig9PV2BQCBs/29tba3y8/N18ODBmLnacaji/T2I999f4j3g94/v31/iPYjk7+84jurq6pSXl6eEhO6PJImKnpOEhAQNGDAgYv9/RkZGXK6QbcX7exDvv7/Ee8DvH9+/v8R7EKnfP5QeExcDYgEAgK8QTgAAgK8kPvDAAw/YLsKmxMRETZs2TUlJUXGEKyLi/T2I999f4j3g94/v31/iPfDb7x8VA2IBAED84LAOAADwFcIJAADwFcIJAADwFcIJAADwFcIJAADwlbgOJ08//bQGDRqk1NRUXXHFFdqwYYPtkjyxePFiXX755UpPT1dWVpZmz56tyspK22VZ88gjjygQCGjBggW2S/HU4cOH9e1vf1v9+/dXWlqaRo0apU2bNtkuyxPNzc267777VFhYqLS0NF1yySV68MEHQ744WTR5++23dcMNNygvL0+BQEAvv/xyu393HEf333+/cnNzlZaWphkzZmj37t2Wqo2Mrt6DpqYm3XPPPRo1apTOO+885eXlac6cOTpy5IjFisPrXOtAWz/4wQ8UCAT0+OOPe1hhq7gNJ7///e+1cOFCLVq0SKWlpRozZoxmzpypY8eO2S4t4tasWaN58+Zp3bp1Ki4uVlNTk6677jrV19fbLs1zGzdu1K9+9SuNHj3adime+vjjjzVlyhT16dNHr732mt5//339/Oc/1wUXXGC7NE88+uijWrJkiZ566int2LFDjz76qH72s5/pF7/4he3SIqa+vl5jxozR008/3eG//+xnP9OTTz6ppUuXav369TrvvPM0c+ZMnTlzxuNKI6er9+DTTz9VaWmp7rvvPpWWlurPf/6zKisrdeONN1qoNDLOtQ64Vq5cqXXr1ikvL8+jyjrgxKmJEyc68+bNa/m5ubnZycvLcxYvXmyxKjuOHTvmSHLWrFljuxRP1dXVOUOGDHGKi4udqVOnOnfeeaftkjxzzz33OFdddZXtMqy5/vrrndtvv73dY9/4xjecW265xVJF3pLkrFy5suXnYDDo5OTkOP/xH//R8tgnn3zipKSkOP/1X/9lo8SI++J70JENGzY4kpz9+/d7VJV3Ovv9Dx065Fx88cVORUWFM3DgQOexxx6zUJ3jxGXPSWNjozZv3qwZM2a0PJaQkKAZM2bovffes1iZHTU1NZKkfv36Wa7EW/PmzdP111/fbj2IF3/5y180YcIE3XTTTcrKytLYsWP17LPP2i7LM5MnT1ZJSYl27dolSdq6davWrl2rWbNmWa7Mjn379qm6urrdZyEzM1NXXHFFXG4TXTU1NQoEAurbt6/tUjwRDAZ166236u6779aIESOs1uKPeWo9duLECTU3Nys7O7vd49nZ2dq5c6elquwIBoNasGCBpkyZopEjR9ouxzMrVqxQaWmpNm7caLsUKz744AMtWbJECxcu1E9+8hNt3LhRP/zhD5WcnKy5c+faLi/ifvzjH6u2tlZDhw5VYmKimpub9dBDD+mWW26xXZoV1dXVktThNtH9t3hz5swZ3XPPPbr55pvj5krFjz76qJKSkvTDH/7QdinxGU7Qat68eaqoqNDatWttl+KZgwcP6s4771RxcbFSU1Ntl2NFMBjUhAkT9PDDD0uSxo4dq4qKCi1dujQuwskf/vAH/e53v9Py5cs1YsQIlZWVacGCBcrLy4uL3x9da2pq0je/+U05jqMlS5bYLscTmzdv1hNPPKHS0lIFAgHb5cTngNgLL7xQiYmJOnr0aLvHjx49qpycHEtVeW/+/Pl69dVXtWrVKg0YMMB2OZ7ZvHmzjh07pnHjxikpKUlJSUlas2aNnnzySSUlJam5udl2iRGXm5ur4cOHt3ts2LBhOnDggKWKvHX33Xfrxz/+sb71rW9p1KhRuvXWW3XXXXdp8eLFtkuzwt3uxfs2UWoNJvv371dxcXHc9Jq88847OnbsmAoKClq2i/v379ePfvQjDRo0yPN64jKcJCcna/z48SopKWl5LBgMqqSkRJMmTbJYmTccx9H8+fO1cuVKvfXWWyosLLRdkqemT5+u8vJylZWVtdwmTJigW265RWVlZUpMTLRdYsRNmTLlrNPHd+3apYEDB1qqyFuffvqpEhLab/4SExMVDAYtVWRXYWGhcnJy2m0Ta2trtX79+rjYJrrcYLJ79269+eab6t+/v+2SPHPrrbdq27Zt7baLeXl5uvvuu/XGG294Xk/cHtZZuHCh5s6dqwkTJmjixIl6/PHHVV9fr+985zu2S4u4efPmafny5XrllVeUnp7eckw5MzNTaWlplquLvPT09LPG15x33nnq379/3Iy7ueuuuzR58mQ9/PDD+uY3v6kNGzbomWee0TPPPGO7NE/ccMMNeuihh1RQUKARI0Zoy5Yt+s///E/dfvvttkuLmFOnTmnPnj0tP+/bt09lZWXq16+fCgoKtGDBAv37v/+7hgwZosLCQt13333Ky8vT7NmzLVYdXl29B7m5ufqHf/gHlZaW6tVXX1Vzc3PLtrFfv35KTk62VXbYnGsd+GIY69Onj3JyclRUVOR1qfF7KrHjOM4vfvELp6CgwElOTnYmTpzorFu3znZJnpDU4e25556zXZo18XYqseM4zn//9387I0eOdFJSUpyhQ4c6zzzzjO2SPFNbW+vceeedTkFBgZOamuoMHjzY+Zd/+RenoaHBdmkRs2rVqg4/93PnznUcx5xOfN999znZ2dlOSkqKM336dKeystJu0WHW1Xuwb9++TreNq1atsl16WJxrHfgim6cSBxwnhqdEBAAAUScux5wAAAD/IpwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABfIZwAAABf+X/94+X6SYjZUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f9d335d3990>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x7f9d33551050>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = c[:m][111];\n",
    "s = c[:s][111];\n",
    "PyPlot.plot(y, linestyle=\"None\", marker=\"+\", color = \"r\")\n",
    "PyPlot.plot(m[s], linestyle=\"-\", marker=\"*\", color = \"b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernel_info": {
   "name": "julia-0.5"
  },
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
